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ABSTRACT 


The  main  problem  of  artificial  satellite  theory  is  a  restricted  two  body  problem 
in  which  the  Legendre  Polynomial  representation  of  the  cylindrically  symmetric 
potential  contains  only  the  first  two  terms.  A  generalized  asymptotic  expansion  is 
used  to  obtain  a  first  order  approximation.  The  solution  at  the  critical  inclination 
is  seen  to  be  of  a  different  type  than  at  other  inclinations.  The  solution  is  finite  for 
all  eccentricities  and  inclinations  when  suitably  restricted  in  time. 
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NOTATION 


q  right  ascension  in  spherical  coordinates 
/?  declination  in  spherical  coordinates 

7  vernal  equinox 
c  cos i o 

e0  eccentricity  of  the  initial  instantaneous  ellipse 
h0  initial  magnitude  of  angular  momentum 
i  inclination 
?o  initial  inclination 

Ji  coefficient  of  second  Legendre  polynomial  in  the  representation  of  the 
gravitational  potential 
J  3/2(JiR2/PS) 

I\t  undetermined  constants  used  in  the  solution 

Po  hl/GM  where  G  is  the  gravitational  constant  and  M  is  the  mass  of  the 
primary  body 

R  equatorial  radius  of  the  primary  body 
r  magnitude  of  the  position  vector  to  the  satellite 
s  sin  ?0 

T  kinetic  energy  of  the  satellite 
t  time 

<o  initial  value  of  time 
u  p0/r 

V  gravitational  potential 

u-’o  initial  value  of  the  argument  of  perigee 

8  angle  measured  from  the  ascending  node  to  the  satellite 
within  the  reference  plane 

0o  initial  value  of  0 

fl  longitude  of  the  ascending  node  measured  from  7 
ft0  initial  value  of  0 
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A.  INTRODUCTION 

The  main  problem  of  artificial  satellite  theory  is  a  restricted  two  body  problem 
in  which  the  Legendre  Polynomial  representation  of  the  cylindrically  symmetric 
potential  contains  only  the  first  two  terms.  Deceptively  simple  in  statement,  the 
main  problem  continues  to  evade  satisfactory  solution.  We  make  no  attempt  to 
survey  the  vast  literature  pertaining  to  the  main  problem,  nor  to  recount  the  history 
of  the  so-called  critical  inclination.  In  fact,  every  effort  is  made  to  make  the  solution 
as  accessible  to  the  non-specialist  as  to  the  astrodynamicist. 

Our  goals  in  attempting  to  obtain  a  solution  are  twofold.  First,  we  desire  a 
solution  employing  a  coordinate  system  based  upon  physical  events  rather  than  an 
abstract  set  of  transformations.  For  instance,  the  use  of  averaged  quantities  such  as 
the  “mean  orbital  plane”  or  other  such  non-physical  artifices  is  eschewed.  Second, 
we  seek  a  solution  which  does  not  tend  to  infinity  at  any  inclination  and  which 
places  no  constraints  on  any  dependent  variable. 

Needed  for  the  solution  are  a  coordinate  system,  asssumptions  concerning 
the  dependent  variables  involved,  and  an  algorithm  for  calculating  the  unknowns. 
These  are  outlined  in  detail  in  the  first  chapters,  after  which  the  solution  is  obtained 
through  extensive  use  of  M  ACSYM  A.  A  discussion  of  the  critical  inclination  problem 
follows. 

Because  the  purpose  and  scope  of  this  work  are  nearly  identical  to  those 
of  Snider  [Ref.  1],  there  are  many  similarities  in  organization,  notation,  and  lan¬ 
guage.  The  principal  difference  in  treatment  lies  in  the  geometry  used  rather  than 
in  approach.  Snider’s  original  approach  is  succesfully  employed  using  an  alternate 
geometry,  and  his  influence  in  this  work  is  pervasive. 
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B.  ORBITAL  KINEMATICS 


Figure  1  shows  (he  reference  system  of  spherical  coordinates  (r.o.  0).  The 
radial  distance  r  is  measured  from  the  center  of  the  planet  O  to  the  satellite  The 
line  0~)  is  in  a  direction  fixed  with  respect  to  an  inertial  coordinate  system.  The 
right  ascension  a  is  the  angle  measured  in  the  planet's  equatorial  plane  eastward 
from  the  line  O 7.  The  declination  or  latitude  0  is  the  angle  measured  northward 
from  the  equator.  The  position  vector  r  of  the  sateiiite  in  the  spherical  coordinate 
system  is 

r  =  r(cos  a  cos  3)bj  +  r(sin  a  cos  5)b2  +  r(sin  J)b3  ( 1 ) 

where  (b1.b2.b3)  are  orthonormal  base  vectors  fixed  in  the  directions  shown. 


polar  axis 


Figure  1:  Spherical  coordinate  system. 

Figure  2  shows  how  we  can  locate  the  satellite  by  its  polar  coordinates  ( r,6 ) 
within  a  rotating  orbital  plane  that  contains  its  position  and  velocity  vectors.  Here 
6  is  the  argument  of  latitude,  i.e.,  the  angle  measured  in  the  orbital  plane  from  the 
ascending  node  to  the  satellite.  The  orbital  plane  is  inclined  at  an  angle  ?  to  the 
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equatorial  plane  and  intersects  the  equatorial  plane  in  the  line  of  nodes,  making  an 
angle  with  the  O')  line. 

orbifol  plone 


Figure  2:  Orbital  plane 

We  introduce  another  orthonormal  set  of  basis  vectors  (B,.Bj.B3)  which 
move  with  the  satellite  so  that  Bj  is  in  the  direction  of  the  position  vector  r.  Bj  is 
also  in  the  orbital  plane  and  B3  =  Bi  x  Bj. 

The  basis  (bj.b2.b3)  may  be  transformed  into  the  basis  (Bj.B2.B3)  bv  a 
succession  of  three  rotations.  First  the  basis  (bj.b^bj)  is  rotated  about  the  b3 
direction  by  the  angle  f).  next  the  basis  is  rotated  about  the  new  first  coordinate 
vector  by  the  angle  ?,  and  finally  the  basis  is  again  rotated  about  the  new  third 
coordinate  vector  by  the  angle  6.  The  two  sets  of  base  vectors  are  related  by  the 
product  of  the  rotation  matrices  representing  each  successive  rotation: 


/  cos  6  sin  0  \  /  1  0  0  \  /  cosfi  sin  fi  0  \  /  b;  \ 

I  -  sin  6  cos 6  0  I  10  cos  i  sin  1  (  -  sin  cosfi  Oil  b?  I  (2) 

\  0  0  1  /  \  0  -  sin  1  cos  1  /  \  0  0  1  /  \  b3  J 
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or 


B, 

B2 

b3 

The  position  vector  r  has  only  one  component  in  the  rotating  basis: 

r  =  rBi  (3) 

Using  the  first  of  Equations  (2),  we  obtain  the  components  of  r  in  the  fixed  basis: 

r  =  r(c os  9  cos  —  sin  9  cos  i  sin  fi)bj 

-f  r(cos  9  sin  +  sin  9  cos  i  cos  Q)b2  +  r(sin  6  sin  ?)b3  (1 ) 

Equating  the  components  of  Equations  (1)  and  (4),  we  can  obtain  the  following 
relations  among  the  angles  (q,^)  of  the  spherical  coordinate  system  and  the  astro¬ 
nomical  angles 


cos  6  cos  ft  -  sin  9  cos  i  sin  Q  cos  6  sin  0.  +  sin  0  cos  i  cos  fi  sin  6  sin  i 
—  sin#  cos  fi  —  cos#  cos « sin  fl  -  sin  <?sin  Q  +  cos  9  cos  icosQ  cost?  sin? 
sin  i  sin  Q  -sintcosQ  cos? 


sin  $  =  sin  9  sin  ? 


(5) 


cos  =  cos  9  sec(a  —  f?) 

The  velocity  dr/di  of  the  satellite  is  obtained  by  differentiating  (3)  with 
respect  to  the  time  t: 

dr  dr  dB] 

Jt=7t*'  +  r~dT  (6) 

Since  the  orbital  plane  must  contain  the  velocity  vector,  we  have  to  enforce 

dBi 


dt 


•  B,  =  0 


(T) 


Substitution  of  Equation  (2)  into  Equation  (7)  leads  to  a  relationship  which  uncou¬ 
ples  the  equations  for  Q {6)  and  i{6): 

dU  tan  9  di 


dO  sin?  dO 


SI 
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The  velocity  (6)  can  then  be  written 


dr  dr  dd  /  .  .  ,n. 

Jt  =  TtB'+rS  +  (9) 

In  the  following  part  of  this  paper,  we  obtain  expressions  for  r(0),  i(0),  ft(0), 
and  dt/d9(9).  The  position  and  velocity  vectors  of  the  satellite  are  then  calculated 
from  the  formulas  in  this  chapter.  The  classical  orbital  elements  p,  e,  and  u  are 
the  semilatus  rectum,  eccentricity,  and  argument  of  perigee  of  the  instantaneous 
(osculating)  ellipse  determined  by  the  position  and  velocity  vectors.  If  needed,  p(0), 
e(0),  and  u>(0)  can  be  obtained  from  our  solution  r(0)  and  dt/d6{9): 


e  cos (6  —  u ;)  = - 1 

r 

esi„(«-w)  =  £(^) 

Numerical  integration  of  the  expression  for  dt/d6{9)  allows  a  transformation  between 
9  and  t.  Figure  3  shows  the  satellite  at  time  t  and  its  relationship  to  the  various 
variables  used. 

This  method  employs  the  true,  rather  than  mean ,  orbital  plane  to  specify  the 
satellite’s  position  and  is  due  to  Struble  [Ref.  2,  3,  4], 

C.  EQUATIONS  OF  MOTION 

The  expressions  in  spherical  coordinates  for  the  kinetic  and  potential  energies 
per  unit  mass  of  a  satellite  orbiting  around  an  oblate  planet  are  respectively: 


r  =  i  [- 

2  \dt 


+  r2  cos2  0  — 


Figure  3:  Satellite  at  time  t. 


where  G  is  the  gravitational  constant,  M  is  the  mass  of  the  planet,  and  R  is  the 
equatorial  radius  of  the  planet.  Substitution  of  these  equations  into  Lagrange’s 
equations 


results  in  the 


i_diT-v)  a 
a  a  (*)  a,(  >- 

following  equations  of  motion: 


q  =  r,  a,  or  l3 


cPr 

It* 


—  r 


d_ 

dt 


2 


dv_ 

dr 


(12) 


d_ 

dt 


+  r2  sin  3  cos  3 
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d\ 

33 


6 


(13) 


Initial  conditions  are  established  by  requiring  that  at  the  initial  time  t0  the 
orbital  parameters  of  the  Keplerian  two  body  ellipse,  determined  by  the  initial  posi¬ 
tion  and  velocity  vectors,  are  known.  The  actual  orbit  is  then  tangent  to  this  initial 
instantaneous  ellipse  (see  Figure  3).  Equating  the  initial  position  and  velocity  vec¬ 
tors  given  by  Equations  (3)  and  (9)  to  the  two  body  expressions,  we  obtain 


r(t0)  = 


Po 


1  +  Co  cos(0Q  —  <^o) 


dr  eoh0s\n(60  -  uQ) 

-_^0)  =  - - - 

at  Po 


(14) 


dd_ 

dt 


(*o) 


ho 

rl  1  +  tan0ocoUo^(0o)] 


(15) 


i{0 o)  =  t0 


(16) 


fl(0o)  =  Oo  (17) 

Here  h0  —  y/GMpo  is  the  initial  value  of  the  satellite’s  specific  angular  momentum 
about  the  center  of  the  planet.  The  subscript  0  on  a  symbol  denotes  that  the 
parameter  is  evaluated  at  the  initial  time  t0. 

We  immediately  have  two  integrals  of  the  equations  of  motion: 

T  +  V  =  constant  (18) 

and 

2  In  do 

r  cos  p—  =  constant  (19) 

dt 

Equation  (18)  states  that  the  mechanical  energy  of  the  satellite  remains  constant. 
Now,  from  Equations  (1)  and  (19) 


dr  ,  22  jdo 

r  x  —  •  b3  =  r  cos  3—  =  h0cosi 0 
dt  dt 


(20) 


i 


Equation  (20)  states  that  component  along  the  polar  axis  of  the  specific  angular 
momentum  of  the  satellite  remains  constant.  Inserting  Equations  (3)  and  (9)  into 
Equation  (20),  and  applying  the  initial  condition  (15),  we  obtain 


dt 

d§ 


r2  cos  i 
h0  cos  zq 


1  +  tan  0  cot  i 


.di\ 
1  dO  ) 


(21) 


This  allows  the  independent  variable  to  be  changed  from  t  to  6. 

Letting  u  =  po/r,  and  using  Equations  (5),  (20),  and  (21),  we  can  rewrite  the 
remaining  equations  of  motion  (12)— (13): 

di  —  2  J a  sin  9  cos  9  sin  i  cos2  i 


d9  -f  2  Ju  sin2  9  cos3  i 

cost 


(22) 


cPu 

lip 


cos2 i  J cos2 i  r 

+  u=  — ^ - 1 - o — 


u2(l  —  3sin2  9  sin2  i) 


du  r.  (Pu  -  , 

+  2u—  sin  0  cos  0(1—3  cos  ?')  —  4 u-j—  sin  9  cos  i  —  2 

do  doJ 


du 

le . 


sin2  9  cos2  i 


(23) 


4  J2u  sin3  6  cos6  i 


u~cos9(2  +  sin2  i)  +  —  (u~  \  sin  6  cos2  i 
dO  d9  \  dO  ) 


The  terms  with  (PujdO 2  on  the  right  side  of  (23)  can  be  eliminated,  yielding  the 
equivalent  equation 


cPu  f 
w  +  u-\ 


2  2  •  t  2  2  • 

C  COS  l  —  J  C  XL  COS  l 


du 

2—  sin  9  cos  0(3  cos2  i  —  1 )  —  u 
du 


—  472u2^-  sin3  0  cos  0  cos6  i( 3  —  cos2  i)  j  (24) 

dO  J 

/(c4  +  4  Juc2  sin2  0cos4  i  +  4u2  J2  sin4  0cos8  i) 

Here  we  have  introduced  the  shorthand  notation  c  =  cos  to.  5  =  sin  ?o<  J  = 

3J2/?72/i 
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The  differential  Equations  (22)-(23)  are  coupled  by  the  nonlinear  terms  and 
apparently  cannot  be  solved  analytically.  If  we  expand  the  right  sides  of  (22)  and 
(24)  in  a  Taylor  series  expansion  in  powers  of  J,  the  equations  simplify  to 
di  —  2Ju  sin  6  cos  0sin  i  cos3  i 


d6 


+  4  J^u^sc3  sin3  0  cos  0  +  0(J3) 


(25) 


d2u 
de 2 


+  U  = 


cos2  i  J  cos2  i  (  —4u  sin2  0  cos4  i 


+ 


-f  u2[l  +  sin2  0(7  cos2  i  —  3)1 


+ 


4  J2u  sin2  6  cos6  i 


3 u  sin2  6  cos4  ? 


du  „  ( du\ 2  o  «  "i 

4-  2u  —  sin  6  cos  0(1—3  cos2  i)  —  2  (  —  )  sin2  6  cos2  i  > 
du  \d6  /  J 

u2[— 1  +  3  sin2  0(1—2  cos2  t)j  + 

w  C~ 

du  ( du  \  ^  \ 

4-u-j^-sin0cos0[4sin2?  +  sin2  0(1  —  3  cos2  i)]  +  (  ~  )  sin20cos2il  +  0(J3) 
du  \  dO  J  J 


(26) 


Using  (8)  together  with  (25)  yields 


—  =  — 2Jucsin2  0  +  4 J2u2c3  sin4  0  +  0(J 3) 
dO 


(27) 


Each  of  the  neglected  terms  in  Equations  (25)-(27)  are  indicated  by  the  O 
symbols.  These  are  terms  which  will  be  multiplied  by  J  to  the  third  power  or  higher. 
Note  that  the  Equations  (25)-(26)  are  identical  to  those  used  as  the  starting  point 
in  the  analysis  of  Eckstein,  Shi,  and  Kevorkian  [Ref.  5]. 

D.  METHOD  OF  SOLUTION 

We  will  use  a  perturbation  technique  to  solve  Equations  ( 25)— ( 27) .  Following 
Erdelyi  [Ref.  6],  we  define  the  order  relations  0  and  o  as  follows.  For  two  real-valued 
functions  f(x)  and  g(x),  we  define  the  relation  /  =  0(g)  if  there  exists  some  real 
constant  £  such  that  |/|  <  £|#j  for  all  x  in  some  domain  of  interest.  Similarly,  we 
define  /  =  o(g)  as  x  — ►  x0  if  for  all  e  >  0,  there  exists  some  neighborhood  of  x0  such 
that  |/|  <  e\g\  within  that  neighborhood. 
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The  series  an(@,  J)Jn  is  said  to  be  a  generalized  asymptotic  approximation 
to  N  terms  of  f(0,J )  with  respect  to  the  scale  {Jn}  as  J  — ►  0  if 

N 

f(6,  J)=Y1  an{0,  J)JN  +  0(JN )  as  J  -»  0 

n= 0 

Let  D  be  some  subset  of  the  real  line  to  be  determined.  We  will  say  that  the 
above  generalized  asymptotic  approximation  holds  uniformly  for  &  £  D  if 

N 

f(0,  J)-J2  J)JN  =  0{JN)  as  J  -  0 

n=0 

uniformly  for  0  6  D.  In  order  to  get  uniformity  it  will  sometimes  be  necessary  to 
bound  the  independent  variable  6 ,  thus  determining  D.  Note  that  these  definitions 
differ  from  the  usual  asymptotic  expansion  definitions  in  that  we  allow  the  coefficient 
functions  a,  to  be  functions  of  both  6  and  J. 

Let  us  suppose  that  each  of  i,  u ,  and  Q  have  formal  generalized  asymptotic 
expansions  for  some  suitably  restricted  values  of  0: 

i*jrtk(9,j)jk  (28) 

k= o 

OG 

u  %  (29) 

*=0 

oo 

fi  *Y,ttk{0,J)jk  (30) 

k=0 

Brenner  and  Latta  [Ref.  7]  found  that  a  modification  was  needed  to  their 
time-like  variable  M  in  order  to  avoid  resonance.  They  accomplished  this  by  means 
of  an  additional  variable  u;  which  multiplied  M,  allowing  secular  term  elimination. 
Analogously  we  introduce  the  variable  y  which  may  be  added  to  integer  multiples  of 
6.  In  this  way  we  too  will  avoid  resonances  except  perhaps  at  certain  inclinations. 
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The  explicit  form  is  then 


y  *52yk{0,J)Jk  (31) 

k—0 

We  further  stipulate  that  for  J  =  0,  u  assumes  the  Keplerian  two  body 
solution  and  that  y  is  the  true  anomaly.  This  forces 

y  ~  +Jyi  +  J2V2  +  ---  (32) 

true  anomaly 

u  «  1  -f  ep  cos  y  +Ju!  +  J2u2  +  •  •  •  (33) 

Keplerian  solution 

An  algorithm  for  the  perturbation  procedure  is  then: 

Let  n—1 

Substitute  Equations  (28),  (30),  (32),  (33)  into  the  equations  of  motion 
Equate  coefficients  of  Jn 
Solve  for  the  nth  order  solution 
Iterate  on  n 

Having  chosen  a  coordinate  system,  made  assumptions  concerning  the  vari¬ 
ables,  and  given  an  algorithm  for  determining  the  unknowns,  we  are  prepared  to 
give  the  solution  in  the  following  chapter. 

E.  SOLUTION 

Substituting  Equations  (33)  and  (28)  into  (25)  and  equating  the  terms  mul¬ 
tiplied  by  J  yields 

=  scsm  20-S-~  sin (y  +  20)  +  ^  sin(y  -  20)  (34) 

A  solution  to  this  equation  is 

*i  =  J  cos  20  H — cos(y  +  20)  +  — ^  cos (y  -  20) 

+  A'i  cos(2y  -  20)  +  I<2  (35) 
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The  last  two  terms  may  be  added  because  they  are,  to  lowest  order,  homogenous 
solutions  to  Equation  (34).  The  term  multiplied  by  the  constant  A'j  was  added  to 
eliminate  secular  terms  in  i2.  Note  that  differentiating  this  term  with  respect  to  0 
produces  terms  multiplied  by  J,  from  Equation  (32).  The  constant  K2  was  added 
to  satisfy  the  initial  condition  (1G),  which  implies  that  i\(9o)  =  0  so 

«*?C  SCCq  SC€-q 

K2  =  —  —  cos  20o - —  cos(30o  —  <^o) - —  cos(0o  +  w0)  ~  Ki  cos  o 

2  6  2 

Substituting  Equations  (28),  (32),  (33),  and  (35)  into  (26),  and  equating 
terms  multiplied  by  J  yields 

<Pu\  3s2  ,  /  5s2  \  l._ 

+  Ul  =  1  _  __  +  +  1  j  +  _[(2  +  5e2)s2  -  2e2]  cos  20 

e2 

+  y(-9s2  +  8)cos2y  +  —(11s2  -  6)cos(y  +  29)  (36) 

4  j 

1^2  r.2  o„ 

+  — 2(3s2-2)cos(2y  +  20)  +  -£(3s2  -  2) - cos(2y  -  20) 

24  8  c 

2sI<2  _L  (odyi  I  A  r  ■ 

+  e°  \2le  + 4  - 5s  j  008 »  +  e° -w sm> 

In  the  above  equation,  the  cosy  and  siny  terms  would  produce  secular  terms  0siny 
and  0cosy  in  ux.  Choosing  dyx/dO  =  5s2 /2  —  2  will  eliminate  these  possibilities. 
Integrating  yields 

2/i  =  - 2  j  (0  -  0O)  +  /\3[ sin(2y  -  20)  +  sin  2u>0]  (37) 

The  term  multiplied  by  K3  was  added  to  eliminate  secular  terms  in  u2.  The  constant 
terms  in  (37)  were  added  to  satisfy  the  initial  condition  y(0o)  =  0o  —  w0. 

A  solution  to  lowest  order  of  Equation  (36)  is  then 

U1  =  1  ~  ~Y  +  e°  (~i“  +  !)  +  +  5eo)  +  2£q] cos 20 

2 

+  pj(9s2  -  8) cos 2y  +  ^(-lls2  +  6)cos(y  +  20) 
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f  t-(— 3s2  +  2)  cos(2 y  +  20)  + 


24 


v(3^2-2)- 


2sKl 


cos(2 y  -  29)  (38) 


2s  I\  i 


+  K4  co s(y  —  29) 


+I< 5  cos (y  -  +  w0)  +  A'6  sin(y  -  90  +  u>0) 


The  term  multiplied  by  K4  was  added  to  eliminate  secular  terms  in  u 7.  The  terms 
multiplied  by  K 5  and  K6  were  added  to  satisfy  the  initial  conditions  (14). 

The  calculations  proceed  by  substituting  Equations  (28),  (32),  (33),  (35), 
(37),  and  (38)  into  (25)  and  equating  terms  multiplied  by  J2: 


di2 

19 


sceg(15s2  —  14) 
24(5s2  -  4)  " 


sin(2y  -  20)  +  . . . 


(39) 


For  brevity  we  have  indicated  on  the  right  side  of  Equation  (39)  only  the  term  that 
would  produce  secular  terms  in  i 2.  Removal  of  this  term  by  making  its  coefficient 
zero  determines  I\\.  Equation  (39)  is  then  integrated  to  determine  i2. 

It  should  be  noted  that  i2  is  needed  for  the  determination  of  the  constants 
I\3  and  K4  in  the  equation  for  u2,  so  integration  is  required. 

Continuing  the  procedure  by  equating  the  terms  multiplied  by  J 2  in  the 
expansion  of  Equation  (26)  determines  y2,  A 3,  and  I\4.  Final  values  of  all  the 
constants  are  listed  in  the  Appendix. 

fl(0)  is  determined  by  substituting  (30)  into  (27)  and  proceeding  in  the  same 
way  as  above.  Note  that  terms  in  J29  arise  in  11(0).  These  must  be  retained  and 
will  restrict  our  variable  0  accordingly. 

In  the  form  below,  use  has  been  made  of  trigonometric  identities  in  order  to 
isolate  terms  containing  the  quantity  (5s2  —  4)  in  the  denominator.  It  may  be  seen 
that  when  (5s2  —  4)  is  zero,  each  of  the  variables  has  a  finite  limit.  Note  the  necessity 
of  keeping  the  J29  terms  in  y  and  fl. 


13 


When  the  quantity  (5s2  —  4)  is  zero,  the  terms  multiplied  by  J  and  J20 
combine,  allowing  the  division  by  (5s2  —  4)  to  occur. 

r  =  p0/|l  +  e0cosj/ +  J  1  -  +  eo  ^1  - 

+  :pr(-(2  +  5eo)s2  +  2ep)cos20  +  ^(9s2  -  8)  cos  2 y 

1  ^  i. 

2 

+  ^(-lls2  +  6)cos(y  +  20)  +  ^(-3s2  +  2)cos(2y  +  20) 
e2 

+  -^(3s2  —  2)  cos(2y  —  20) 

8 

ec[15(2  -(-  e2)s4  -  14(4  +  e§)s2  -4-  24]  sin  \j0  (~  -  2)1  sin[0  +  u;0] 

+ - — - 

e2s2(15s2  -  14)  sin  [j8  —  2)]  sin  [2u0  -  J8  -  2)] 

6(5s2  -  4) 


+ 


2  2  2 

—  cos (y  —  0$  -f  3u.’o)  +  ^t(3s2  —  2)  cos(y  —  3 60  +  3^’o) 
lb  24 


els2 


— -  cos (y  -  560  +  3cj0)  +  ^(3s2  -  2)  cos(y  -  2 0o  +  2u ;0) 
lb  4 


3e0s 


eot 


g  cos(y  -  40o  +  2u;0)  -  +  1) cos(y  +  2u0) 


(40) 


1 


+  ^[(-2  +  5eo)s2  -2ejj]cos(y  +  0O  +  u0) 

+  ^  [6  +  5eg)s2  —  4(1  4-  to)]  cos(y  —  0O  +  u-’o) 

+  ^[-(14  +  5e2)s2  +  2e2]  cos (y  -  30o  +  u>0) 

2  2 

+  ~^(9s2  —  4)  cos  (y  +  30o  —  w’o)  -{-  -^(— 7s2  +  6)  cos(y  +  0o  —  <-^o) 
4o  o 

+  Tr(-5s2  +  4)  cos(y  -  0O  -  w0)  +  ^(2 s2  -  1 )  cos (y  +  20o) 
lb  4 

+  ^(-3s2  +  l)cos(y  -  20o)  +  ^(-3s2  +  2) cosy 
4  4 

eos2 

+  eos2  cos(0o  +  u?o)  +  — ~ —  cos(30o  ~  ^o) 


+  s2  cos  20o 


+  O(J2,J30,...)} 
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y  =  0-u,o  +  ./{(^-2)(0-0o) 

e2(-75s6  +  260s4  -  296s2  +  112)  sin  [j8  -  2)]  cos  [2o,0  -  JO  -  2)J 

+  24(5s2  -  4)2 

,  r2flf€oS2(-15s2  +  14)(15s2  -  13)cos2o,0  e0s2,lc  2  lox  __ _/a  ,  , 

+  J0[ - 24(5 ~s^~T) - +  ~ (155  ~  13)cos(0o  +  wo) 

2  2 

+  yy(15s2  —  13)cos(30o  —  o,0)  +  y(15s2  —  13)cos20o  (41) 


+ 


^[5(9e2  +  34)s4  +  4(9e2  -  34)s2  -  56e2]|  +  0{J\  J36, . . .) 


i  =  io  +  scJ |  ^  cos  20  +  y  cos{y  +  20)  +  y  cos(y  —  20) 

e2(-15s2  +  14)  sin  \j6  -  2)]  sin  [2o,0  -  J6  -  2)] 

+  12(5s2  —  4) 

~  2  cos“^°  — g”  cos(30q  —  w0)  — —  cos(0o  +  o,q)|  +  0(J 2,  J30, . . .) 


(42) 


U  =  f20  +  cJ^0o  —  6  +  ^  sin  20  —  e0sin  y 

+  T"  sin(j/  +  20)  —  —  sm(y  —  20)  —  ^  sin  20o 
6  2  2 

+  e0  sin(0o  -  o,0)  -  ^sin(30o  -  o,0)  -  ^  sin(0o  +  o,0) 

6  2 

e2(15s4  -  45s2  +  28)  sin  [j0  -  2)j  cos  [2u>0  -  JO  -  2)] 

+  6(5s2  -  4)2 

+  cJ29^  0lo((55g2-l4)4)  cos 2u,0  -  e0s2  cos(0o  +  u>0)  (43) 

2  2 
-  yy  cos(30o  -  u>o)  -  S2  cos  20o  +  ^|(7s2  ~  4) 

+  ^(-sJ  +  6)}  +  0(7!,J39,...) 
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t  =  to  + 


vJA 


1  +  J 


(~3s2  +  2) 


cos  20 


2  ,  \ _  .  eo(-4s2  +  3) 


+  e0(s  -l)cosy  + 
eo(—2s2  +  1 ) 


cos(y  -f  20) 


+ 


+ 


cos(y  —  20) 


(44) 


eQS2(1552  -  14)  sin  ^ JO  -  2)]  sin  2*jq  —  JO  —  2) 


12(5s2  —  4) 


e05 


e0s‘ 


+  s  —  1  +  —  cos  20o  +  ~ —  cos (30o  —  wo)  H — r —  cos(0o  +  u?o) 


+  0(J2,  J30, . . .)  j 


do 


To  check  the  solution,  we  can  see  if  the  specific  mechanical  energy  (18)  of 
the  satellite  remains  constant.  Substitution  of  the  solution  (40)-(42)  into  Equations 
(10)  and  (11 )  yields 

T  y  G-W(l-eg)  G.W,fl*(l-3smH)  ,  , 

~Po  2|r(i0)P 

The  first  two  terms  on  the  right  side  are  easily  recognized  as  the  value  of  the  specific- 
mechanical  energy  at  the  initial  time  tQ. 

As  a  further  check  on  the  solution,  we  can  see  if  it  reduces  to  previous  results 
for  equatorial  and  polar  orbits,  obtained  by  Danielson  and  Snider  [Ref.  8].  Setting 
t’o  =  0  and  using  the  independent  variable  a  measured  from  the  line  0~),  we  find  that 
Equations  (4 0)— (42)  reduce  to  Equations  ( 1 8)— (22)  in  [Ref.  8].  Setting  i0  =  tt/2  and 
using  the  expansion  cos(y  +  Jk)  tss  cos  y  —  Jk  sin  y,  we  find  that  Equations  (40)— (42) 
reduce  to  Equations  ( 38 )— (4 1 )  in  [Ref.  8].1 

'There  are  misprints  in  Equations  (37)— (40)  of  [8].  The  term  2e0c2  cos  y  should  be  added  on  the 
right  of  Equation  (37).  The  third  from  last  term  in  Equation  (38)  should  be  eo/3cos(3Jo-~o)  The 
term  — c272J  should  be  added  on  the  right  of  Equation  (39).  The  sign  of  the  second  trigonometric 
term  in  the  expression  for  c$  should  be  changed:  (— e0/12  +  eo/32)sin(  Jo  +  -'o) 
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F.  THE  CRITICAL  INCLINATION 


There  exist  many  “solutions”  to  the  main  problem.  Most  of  these  are  un¬ 
satisfactory.  Some  are  so  abstruse  as  to  have  little  or  no  physical  meaning,  while 
others  suffer  from  having  a  restricted  domain  in  inclination,  eccentricity,  or  both. 
Chief  among  the  difficulties  with  proposed  solutions  is  their  behavior  at  certain  in¬ 
clinations.  The  most  persistent  of  these  troublesome  inclinations  is  the  so-called 
critical  inclination  icr *t  =  ±arcsin2/v/5  (mod  tt).  This  manifests  itself  in  the  form 
of  denominator  terms  like  (4  —  5sin2z),  (1  —  5cos2i),  or  (tan  i  —  2).  When  2  as¬ 
sumes  this  critical  inclination,  2cm ,  each  of  these  terms  become  zero,  rendering  the 
solutions  useless.  It  is  at  this  point  that  we  must  confront  the  divisor  (5s2  —  4)  in 
Equations  (40)— (44 ).  As  remarked,  each  of  (40)-(44)  has  a  finite  limit  at  the  critical 
inclination.  How  does  this  come  about?  At  the  critical  inclination,  the  quotient 
sin(J0(5s2/2  —  2)/(5 s2  -  4)  is  replaced  by  the  limit,  yielding  a  term  in  JO.  This,  in 
the  parlance  of  orbital  mechanics,  is  an  odious  secular  term  —  an  unbounded  term 
which  grows  with  the  time-like  variable.  The  expressions  (40)-(44)  at  the  critical 
inclination  are: 


f  3s2  (  5s2  \ 

=  p0/ 1 1  +  e0  cos  y  +  J  1  -  —  +  e2  I  1  -  —  I 


+  +y^(-(2  +  5c2)s2  +  2c2)  cos  20  +  ^(9s2  -  S)  cos2.y 


+  ^(  — lls2  +  6)cos(y  +  20)  +  ^(— 3s2  +  2)cos(2y  +  20) 

24  24 

+  |(3s2-2)cos(2y-20) 

2  2  2 

—  cos(y  —  0o  +  3u;o)  +  -^-(3s2  —  2) cos (y  —  3 0O  +  3^'o) 
10  24 


els2 


—  -  cos (y  —  50o  +  3w'o)  +  —(3s2  —  2)  cos (y  —  20o  +  2^’o) 

16  4 


3e0s 


fo , 


cos {y  -  400  +  2u,’o )  -  ^(s2  +  1 )  cos (y  +  2~-0) 
o  4 


(45) 
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+  ^[(—2  +  5eo)s2  -  2eo]cos(y  +  do  +  ^o) 

+  ^[(6  +  5eo)s2  -  4(1  +  eo)]cos(y  -  0O  +  ^o) 

+  rr[— (14  +  5cq )s2  +  2c\)  cos (y  -  30o  +  ^o) 

£t*T. 

e1  e2 

+  t|(952  -  4)  cos (y  +  3 60  -  w0)  +  -£(~7s2  +  6)  co s(y  +  0o  -  u,’0) 

48  o 

2 

+  ^-(-5s2  +  4)  cos(y  -  0o  -  u>0)  4  y(2s2  ~  1 )  cos(y  4  2 0O) 
lb  4 

+  —(—3s2  4  1)  cos(y  -  20o)  4  -^(-3s2  4  2)  cos  y 
4  4 

e2s2 

-f  e0s2  cos(0o  +  ^o)  +  ~^r~  cos(30o  —  <^o)  4  s2  cos  2 0o 

o 

4  J20f^(15(2  +  e2)s4  -  14(4  4e2)s2  4  24)sin(0  4u,o) 

.24 

4  ^(15s2-  14)sin(2a,0)]  +  0( J2,  J30, . . .)} 

y  =  0  —  u)o  4  J  |  - 2^  (6  —  #o) | 

2  2 

4  J2^(^(-105s4  4  130s2  -2S)cos2w04  ~(15s2  -  13)  cos(0o  4  ->o) 

1 48  - 

4  —  (15s2  -  13)cos(30o-wo)+  ^(15s2  -  13)cos20o  (‘16) 

6  2 

4  —  [5(9e2  +  34)s4  +  4(9e2  -34)s2-56e2]  )  4  0(J2,  J30, . . .) 

96  L  > 

i  -  z0  4  sc,/ cos  20  +  ^cos(y  4  26)  4  -^cos(y  -  20) 

-  -  cos  20o  -  -§■  cos(30o  -  <^o)  -  ir  cos(0o  4  ^o) 

2  b  2 

4-  J2^f££2(_i5s2  -f  14)  sin  2u>o  4  0{J2,  J30, . . .) 


IS 


(48) 


fi  =  fi0  +  cjjflo  —  $  +  ^  sin  20  —  e0sin  y 

+  sin(j/  +  20)  —  ^  s\n(y  —  20)  —  ^  sin  20o 
b  l  l 

+  e0  sin(0o  -  w0)  —  ^  sin(30  -  w0)  -  ^  sin(0o  —  w0) 
b  2 

(  e2 

+  cJ20l-~(6s2  —  7)cos(2u0)  —  eo'S2cos(0o  4-  w0) 

2  2 
-  -77-  cos(30o  -  w0)  -  s2  cos  20o  +  p(7s2  -  4) 

O 


t  =  t 


„4/V{ 

n o  JBq  l 


1  +  J 


(—3s2 +  2) 


cos  20 


2  „ _  ,  e0(-4s2  +  3) 


+  eo^  -  l)cosy  + 
e0(— 2s2  +  1) 


cos(y  +  20) 


cos  {y  —  29) 

/ 

^2  ^  ^2  ^0*5^ 

+  s2  —  1  +  —  cos  20o  H — — —  cos (30o  —  Wo)  H - - —  cos(0o  +  u.’o) 

lb  2 


(49) 


2  2 

+  J26^-{\bs2  -  14)  sin(2u.’0)  +  0(J2,  J30, . . .)} 
24  J 


dd 


As  can  be  seen,  none  of  the  variables  fails  to  be  defined  at  any  inclination  or 
eccentricity. 

A  long-standing  debate  in  astrodynamics  centers  about  the  nature  of  the 
critical  inclination’s  ubiquitous  presence  in  proffered  solutions.  Of  course  something 
important  does  occur  at  the  critical  inclination  —  the  line  of  apsides  remains  fixed. 
That  this  event  should  also  cause  most  attempts  at  solving  the  problem  to  fail  has 
been  the  subject  of  much  debate  and,  apparently,  misunderstanding. 

Examination  of  Equation  (45)  shows  that  when  i  assumes  the  critical  incli¬ 
nation  there  result  terms  in  J20.  At  the  order  to  which  we  have  approximated  the 
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solution  the  form  of  our  approximation  at  the  critical  inclination  differs  from  that  of 
all  inclinations.  In  fact,  if  i  ^  it  appears  that  we  can  obtain  arbitrary  accuracy 
for  sufficiently  small  J  and  6  in  some  J  dependent  interval  [0,c/Jm],  (m  constant, 
integer,  c  an  0(1)  constant),  by  taking  sufficiently  many  terms.  However,  when 
*  =  fcrit,  the  expansion  has  a  restricted  range  of  validity,  for  6  in  [0,  c/J]  where  c  is 
an  0(1)  constant.  In  other  words,  while  we  have  dropped  the  higher  order  terms  and 
higher  order  secular  terms  in  our  solution,  it  is  hoped  that  away  from  the  critical 
inclination  it  will  be  possible  to  improve  the  approximation  to  arbitrary  accuracy 
simply  by  continuing  the  method  used.  Not  so  at  the  critical  inclination — we  are  al¬ 
ready  confronted  with  a  secular  term  which  cannot  be  eliminated,  regardless  of  how 
many  terms  we  keep  in  our  approximation.  We  are  thus  in  agreement  with  Coffey, 
Deprit,  and  Miller  [Ref.  9]  that  the  critical  inclination  is  an  intrinsic  singularity, 
independent  of  the  method  used  to  solve  the  problem.  Moreover,  there  may  emerge 
other  critical  inclinations  as  higher  order  analyses  are  performed.  It  may  be  that 
a  higher  order  approximation  will  yield  other  inclinations  where  secular  terms  arise 
which  cannot  be  eliminated.  The  well-known  critical  inclination  might  be  only  the 
first  in  a  sequence. 

Enters  now  the  astrodynamicist,  the  man  in  the  field.  When  the  equations 
of  motion  are  integrated  at  the  critical  inclination,  no  strange  effects  are  reported. 
The  mathematicians  have  declared  the  critical  inclination  a  glaring  singularity,  yet 
everyday  practitioners  consider  it  a  trifling  locale.  If  a  satellite  were  to  start  an 
orbit  at  the  critical  inclination,  Equation  (45)  predicts  a  secular  variation  in  u(6) 
of  order  J26.  Clearly  such  slow  growth  with  6  explains  why  no  appreciable  effect  of 
the  critical  inclination  has  been  encountered  in  the  intensive  numerical  integrations 
which  have  been  carried  out  over  the  years.  Additionally,  if  a  satellite  is  initially  at 
the  critical  inclination,  Equation  (47)  states  that  the  expression  for  its  inclination 
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also  contains  a  term  in  J26.  This  means  that  the  satellite  will  tend  away  from  the 
critical  before  the  secular  terms  which  arise  in  u(0)  at  the  critical  become  important. 
As  the  satellite  moves  away  from  the  critical  inclination,  the  terms  in  u{9)  resume 
their  bounded  trigonometric  form.  From  a  practical  viewpoint,  the  critical  inclina¬ 
tion  problem  is  in  fact  a  nonproblem.  The  simplifications  and  assumptions  made 
will  perforce  limit  the  time  over  which  the  solution  is  valid,  further  obscuring  at¬ 
tempts  to  physically  observe  anomalies  at  the  critical  inclination.  We  must  concede 
that  other  perturbative  forces  would  play  a  major  role  in  determining  the  degree  to 
which  the  critical  inclination  could  induce  observable  effects.  This  is  not  to  endorse 
a  variety  of  patchwork  measures  now  being  used  to  determine  satellite  motion,  such 
as  “gapping”  the  critical  inclination  in  order  to  avoid  dividing  by  zero.  The  critical 
inclination,  delicate  as  though  it  may  be,  is  perhaps  a  powerful  telltale  for  faulty 
modeling. 

G.  CONCLUSIONS 

As  shown,  our  solution  satisfies  the  goals  given  in  the  introduction. 

If  the  perturbation  procedure  is  continued,  it  is  anticipated  that  each  of  the 
additional  terms  added  to  Equations  (40)-(44)  will  be  multiplied  by  one  of  the 
factors  (J2,  J3,  J30,  J40, . . .).  If  we  restrict  0  <  1  /«/,  the  neglected  terms  should  be 
order  J 2.  (For  an  Earth  satellite  J  <  3/2  x  1 0-3,  so  for  at  least  100  revolutions  the 
relative  error  should  be  order  10~6.) 

If  we  restrict  6  <  1,  all  of  the  terms  of  order  J20  in  the  solution  (40)-(44)  can 
be  dropped  without  increasing  the  order  of  magnitude  of  the  error.  This  considerably 
simplifies  the  solution  to: 
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APPENDIX 
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